clear; clc; close all;

%% 固定参数
a0 = -2;
a1 = 0.83;
a2 = 0.14;
b2 = 1;
d1 = 0;
d2 = 0;

% 扫描参数 b0 和 b1 的取值范围
b0_values = linspace(-5, 5, 101);   % 例如：-3 到 3，共31个点
b1_values = linspace(-1, 1, 101);     % 例如：-1 到 1，共31个点
nb0 = length(b0_values);
nb1 = length(b1_values);

% 迭代设置
N_total = 1000;    % 总迭代步数
N_trans = 200;     % 舍弃瞬态步数
N_show  = N_total - N_trans;  % 有效记录步数

% 初始条件
x0 = 0.1;
y0 = 0.1;
z0 = 0.1;

%% 预分配存储谱熵矩阵, 尺寸为 (nb1, nb0)
% 约定外层循环为 b0, 内层为 b1
SE = zeros(nb0, nb1);

%% 双重循环，遍历每个 (b0, b1) 组合
for i = 1:nb1
    for j = 1:nb0
        b1 = b1_values(i);
        b0 = b0_values(j);
        
        % 初始化状态
        x = x0; y = y0; z = z0;
        
        % 舍弃瞬态
        for n = 1:N_trans
            [x, y, z] = vl_dmsin(a0, a1, a2, b0, b1, b2, d1, d2, x, y, z);
        end
        
        % 记录后续 N_show 步的 x 值到 data 数组
        data = zeros(1, N_show);
        for n = 1:N_show
            [x, y, z] = vl_dmsin(a0, a1, a2, b0, b1, b2, d1, d2, x, y, z);
            data(n) = x;
        end
        
        % 计算归一化谱熵，并存入 SE(i,j)
        % 使用 SpectralEntropy 函数，改进后归一化到 [0,1]
        SE(i,j) = SpectralEntropy(data);
    end
    fprintf('%.1f%%\n', 100 * i / nb1);
end

%% 生成绘图所需的网格
[B1, B0] = meshgrid(b1_values, b0_values);

%% 绘制 Fig.19(a)：归一化谱熵 3D 曲面图
figure;
surf(B1, B0, SE);  % b0 作为 x 轴，b1 作为 y 轴，谱熵作为 z 值
shading interp;
xlabel('b_0','FontSize',12);
ylabel('b_1','FontSize',12);
zlabel('归一化谱熵','FontSize',12);
title('归一化谱熵, b_0, b_1','FontSize',14);
colorbar;


%% 辅助函数： SpectralEntropy
function se_norm = SpectralEntropy(data)
    % 计算归一化谱熵，归一化到 [0,1]
    % data: 输入时间序列 (行向量)
    % se_norm: 归一化后的谱熵值
    %
    % 计算过程：
    %   1. 对数据去均值，然后做 FFT 得到能量谱
    %   2. 只取正频率部分，归一化为概率分布 P
    %   3. 计算谱熵 SE = -sum(P * log2(P))
    %   4. 归一化：SE_norm = SE / log2(M)
    
    data = data(:)';  % 转为行向量
    N = length(data);
    
    % FFT 计算（去均值后）
    Y = fft(data - mean(data));
    P2 = abs(Y).^2;
    % 取前一半频率点
    P1 = P2(1:floor(N/2));
    % 归一化为概率分布
    P_norm = P1 / sum(P1);
    
    % 计算谱熵（采用 log2）
    idx = P_norm > 1e-12;
    SE = -sum(P_norm(idx) .* log2(P_norm(idx)));
    
    % 最大可能熵为 log2(M)，其中 M 为有效频率点个数
    M = length(P1);
    se_norm = SE / log2(M);
end
